home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / dft.stu < prev    next >
Text File  |  1995-03-23  |  9KB  |  273 lines
  1. Article 2502 of comp.sys.handhelds:
  2. Path: en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!uakari.primate.wisc.edu!sdd.hp.com!usc!apple!agate!shelby!portia.stanford.edu!elaine3.stanford.edu!mcgrant
  3. From: mcgrant@elaine3.stanford.edu (Michael Grant)
  4. Newsgroups: comp.sys.handhelds
  5. Subject: A convenient little set of DFT functions, PART 1
  6. Summary: My homework is one heck of a lot easier now
  7. Keywords: DFT FFT DHT FHT CIA NSA ASPCA
  8. Message-ID: <1990Nov21.191129.5321@portia.Stanford.EDU>
  9. Date: 21 Nov 90 19:11:29 GMT
  10. Sender: mcgrant@portia.stanford.edu
  11. Followup-To: comp.sys.handhelds
  12. Organization: Stanford University - AIR
  13. Lines: 256
  14.  
  15. Many of you asked for my programs concerning 1) matrix resolvents,
  16. and 2) DFT's and convolutions of discrete functions.  The former
  17. shall come at a future time, but, for now, here are the DFT 
  18. programs.
  19.  
  20. All of these functions operate on column vectors.  Originally,
  21. they worked on lists, but there are so many functions that 
  22. can be performed on the elements of an array that I just HAD
  23. to convert them.  SO, they don't follow conventional notation
  24. like they used to, but they are now part of a much more
  25. powerful repertoire of functions.
  26.  
  27. The following programs are a set that I have compiled for my
  28. OWN benefit to perform some simple DFT analysis for academic
  29. functions.  They are by no means awesome; they only handle one
  30. dimension, and they aren't terribly fast.  BUT, they sure have
  31. saved time on homework and exams...
  32.  
  33. While I THINK they are bug-free, I had to type them in by hand,
  34. so who knows.  Of course, I don't musspelll tu muche when I type,
  35. so perhaps these will be OK.  In any case, I would appreciate it
  36. if you notify me of any corrections that need to be made
  37. (at mcgrant@portia.stanford.edu).
  38.  
  39. NOTE:  I will post an un-commented set of these next, so that
  40. you will have less editing to do.  Note that I have a 28S, so
  41. these programs may have to be tweaked...
  42.  
  43. ---
  44. I recommend this order in your directory:
  45. { DFT IDFT DHT IDHT CCNV LCNV PROD MP->C C->MP
  46.   PAD2 RPAD ZPAD SAMP AR01 AR02 PRCSN } ORDER
  47.  
  48. 'DFT' -- performs a Discrete Fourier Transform
  49.          on the given input vector.
  50. INPUT: 1: a column vector L
  51. OUTPUT: 1: DFT(L)
  52. REQUIRES: DHT
  53. COMMENTS: This program uses the Discrete Hartley Transform.
  54.           The DHT is a lot faster than the DFT, and the 
  55.           time required to extract the DFT from the DHT is
  56.           minimal in comparison.  Of course, for real speed,
  57.           perhaps someday I will write an FHT/FFT!
  58. << DHT DUP ARRY-> SZ
  59.    << 2 SZ LIST-> DROP 1 - 
  60.       FOR I I ROLL
  61.       NEXT SZ ->ARRY
  62.    >> DUP2 - (0,1) * - +
  63.    IF DUP IM CNRM NOT
  64.    THEN RE
  65.    END .5 *
  66. >>
  67.  
  68. 'IDFT' -- performs an Inverse Discrete Fourier Transform
  69. INPUT: 1: a column vector L
  70. OUTPUT: 1: IDFT(L)
  71. REQUIRES: DFT
  72. COMMENTS: Nothing more than INV(N)*CONJ(DFT(CONJ(L)).
  73. << CONJ DFT CONJ DUP SIZE LIST-> DROP / >>
  74.  
  75. 'DHT' -- performs a Discrete Hartley Transform
  76. INPUT: 1: a column vector L
  77. OUTPUT: 1: DHT(L)
  78. REQUIRES: PRCSN
  79. COMMENTS: The Discrete Hartley Transform uses cas(theta)=
  80.           cos(theta)+sin(theta) as its kernel.  Some people
  81.           will notice that this can be calculated as sqrt(2)*
  82.           sin(theta+pi/4), but for some reason this calculation
  83.           just wasn't quite as accurate as sin+cos; besides,
  84.           the speed savings is not all that great.
  85.           The PRCSN variable tells the DHT to round the answers
  86.           to PRCSN digits.  I do this because 4.0000000001 takes
  87.           up the whole screen, tho' I know it's 4!
  88.           Notice that ~pi~ is the pi key .
  89. << DUP SIZE LIST-> DROP RCLF -> L N F
  90.    << RAD PRCSN FIX L ~pi~ ~pi~ + -> S
  91.       << 0 N 1 - 
  92.          FOR J J N / S * DUP SIN SWAP COS + J 1 + SWAP PUT
  93.          NEXT
  94.       >> -> C
  95.       << L 0 N 1 -
  96.          FOR J 0 0 N 1 -
  97.             FOR K J * N MOD 1 + GET 'L' K 1 + GET * +
  98.             NEXT RND J 1 + SWAP PUT
  99.          NEXT
  100.       >> F STOF
  101.    >>
  102. >>
  103.  
  104. 'IDHT' -- performs an Inverse Discrete Hartley Transform.
  105. INPUT: 1: a column vector L
  106. OUTPUT: 1: IDHT(L)
  107. REQUIRES: DHT
  108. COMMENTS: DHT is symmetric, so this is easy (N^-1*DHT(L)).
  109. << DHT DUP SIZE LIST-> DROP / >>
  110.  
  111. 'SAMP' -- generates an array by sampling a function.
  112. INPUT: 2: a self-contained function F
  113.        1: the number of samples to take
  114. OUTPUT: [ F(0) F(1) F(2) ... F(N-1) ]
  115. COMMENTS: Sorry this isn't more complex guys and gals, but
  116.           it works fine for me.
  117. << -> F N
  118.    << 0 N 1 - 
  119.       FOR J J F EVAL
  120.       NEXT { N } ->ARRY
  121.    >>
  122. >>
  123.  
  124. 'ARO2' -- performs a binary operation on the elements of
  125.           two vectors (thereby combining them).
  126. INPUT: 3: column vector L1 (size S1)
  127.        2: column vector L2 (size S2)
  128.        1: a self-contained binary function F
  129. OUTPUT: [ F(L1(1),L2(1)) F(L1(2),L2(2)) ...
  130.           F(L1(min(S1,S2)),L2(min(S1,S2))) ]
  131. COMMENTS: Notice that, if the two lists are of different
  132.           sizes, that AR02 truncates (in effect) the
  133.           longer one to match the size of the shorter one.
  134. << -> L1 L2 PR
  135.    << L1 SIZE LIST-> DROP L2 SIZE LIST-> DROP MIN -> N
  136.       << L1 { N } RDM 1 N
  137.          FOR J 'L1' J GET 'L2' J GET PR EVAL J SWAP PUT
  138.          NEXT
  139.       >>
  140.    >>
  141. >>
  142.  
  143. 'PROD' -- Performs an element-by-element multiplication.
  144. INPUT: 2: column vector L1 (size S1)
  145.        1: column vector L2 (size S2)
  146. OUTPUT: [ L1(1)*L2(1) L1(2)*L2(2) ... L1(min(S1,S2))*L2(min(S1,S2)) ]
  147. REQUIRES: ARO2
  148. COMMENTS: See ARO2.  This function is useful because convolution
  149.           in one domain is equivalent to multiplication in the
  150.           other domain...but this is not QUITE true for Hartley...
  151. << << * >> ARO2 >>
  152.  
  153. 'AR01' -- Performs a function on each element of the vector.
  154. INPUT: 2: column vector L (size N)
  155.        1: self-contained function F
  156. OUTPUT: 1: [ F(1) F(2) ... F(N) ] 
  157. COMMENTS: well, I just use this for C->MP and MP->C below.
  158. << -> L PR
  159.    << L SIZE LIST-> DROP -> N
  160.       << L 1 N
  161.          FOR J 'L' J GET PR EVAL J SWAP PUT
  162.          NEXT
  163.       >>
  164.    >>
  165. >>
  166.  
  167. 'C->MP' -- converts a vector from complex to mag/phase.
  168. INPUT: 1: column vector L (size N)
  169. OUTPUT: 1: [ R->P(L(1)) R->P(L(2)) ... R->P(L(N)) ]
  170. REQUIRES: AR01
  171. COMMENTS: Convenient when you just want the magnitude
  172.           of your DFT (just do a C->MP RE)!
  173. << (1,0) * << R->P >> ARO1 >>
  174.  
  175. 'MP->C' -- converts a vector from mag/phase to complex.
  176. INPUT: 1: column vector L (size N)
  177. OUTPUT: 1: [ P->R(L(1)) P->R(L(2)) ... P->R(L(N)) ]
  178. REQUIRES: AR01
  179. COMMENTS: When you're too short-sighted to remember to
  180.           save the original answer to the DFT...
  181. << (1,0) * << P->R >> ARO1 >>
  182.  
  183. 'CCNV' -- performs the circular convolution of two vectors
  184. INPUT: 2: column vector L1 (size N)
  185.        1: column vector L2 (size N)
  186. OUTPUT: 1: L1 C* L2 (size N)
  187. COMMENTS: This program will generate an error if the two lists
  188.           are of different sizes.  So, use ZPAD, PAD2, or
  189.           RPAD to adjust the sizes of your vectors.
  190. << DUP SIZE LIST-> DROP -> L1 L2 N
  191.    << IF N L1 SIZE LIST-> DROP <>
  192.       THEN [ 1 ] TRN
  193.       ELSE L1 1 N
  194.          FOR J 0 1 N
  195.             FOR K 'L1' K GET 'L2' J K - N MOD 1 + GET * +
  196.             NEXT J SWAP PUT
  197.          NEXT
  198.       END
  199.    >>
  200. >>
  201.  
  202. 'LCNV' -- performs the linear convolution of two vectors
  203. INPUT: 2: column vector L1 (size N1)
  204.        1: column vector L2 (size N2)
  205. OUTPUT: 1: L1 * L2 (size N1+N2-1)
  206. REQUIRES: PAD2 CCNV
  207. COMMENTS: Gee, the REQUIRES line is longer than the program!
  208. << PAD2 CCNV >>
  209.  
  210. 'PAD2' -- pads two vectors so L1 C* L2 <==> L1 * L2
  211. INPUT: 2: column vector L1 (size N1)
  212.        1: column vector L2 (size N2)
  213. OUTPUT: 2: column vector L1 (size N1+N2-1)
  214.         1: column vector L2 (size N1+N2-1)
  215. COMMENTS: This program is essential to the proper operation of
  216.           LCNV, because circular convolution must 'see' enough
  217.           padded zeros so that no circular overlap occurs.
  218. << DUP2 SIZE LIST-> DROP SWAP SIZE LIST-> DROP + 1 - 1 ->LIST -> S
  219.    << S RDM SWAP S RDM SWAP >>
  220. >>
  221.  
  222. 'ZPAD' -- adds as many zeros as is necessary to extend the
  223.           input vector to at least N elements
  224. INPUT: 2: column vector L (size M)
  225.        1: the desired new vector size N
  226. OUTPUT: IF N>M, 1: RDM(L,N)
  227.         IF N<M, 2: L
  228. COMMENTS: Not much to this.
  229. << -> L N
  230.    << L 
  231.       IF L SIZE LIST-> DROP N <
  232.       THEN { N } RDM
  233.       END
  234.    >>
  235. >>
  236.  
  237. 'RPAD' -- replicate the vector to hit N vectors.
  238. INPUT: 2: column vector L (size M)
  239.        1: the desired new vector size N
  240. OUTPUT: IF N>M: [ L(1) L(2) ... L(M) L(1) .. L(M) ... ] (size N)
  241.         IF N<M: L
  242. COMMENTS: This is useful when circularly convolving two functions
  243.           one's period is a multiple of the other's.  Otherwise,
  244.           I can't think of any mathematically safe use for this.
  245. << -> L N
  246.    << L SIZE LIST-> DROP -> M
  247.       << L
  248.          IF M N <
  249.          THEN { N } RDM M 1 + N
  250.             FOR J 'L' J M MOD
  251.                IF DUP NOT
  252.                THEN DROP M
  253.                END GET J SWAP PUT
  254.             NEXT
  255.          END
  256.       >>
  257.    >>
  258. >>
  259.  
  260. 'PRCSN' -- # of digits to RND to in DHT
  261. COMMENTS: I recommend a fairly large number, though my
  262.           experience is that the DHT is accurate to about
  263.           10 digits (as far as roundoff error goes).
  264. 10
  265.  
  266.  
  267. ---
  268. Michael C. Grant              "God does not play dice." Einstein
  269. Information Systems Lab       "Geez, He'd win a lot if he did,
  270. mcgrant@portia.stanford.edu    though..." Mike
  271.  
  272.  
  273.